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I,  OVERVIEW  AND  SUMMARY 


Fundamental  to  this  work  is  the  development  of  a  continuum 
formulation  that  can  accurately  account  for  the  effects  of 
interlaminar  shear  and  interlaminar  normal  stress  variation 
thru-the-thickness  of  a  laminate.  Furthermore,  emphasis  is 
particularly  on  tapered-twisted  airfoil  geometries  which  can  be 
analytically  represented  as  an  assemblage  of  thin  to  moderately  thick 
finite  elements.  To  achieve  solution  efficiencies,  the  elements 
developed  in  this  work  are  of  the  triangular/quadrilateral  plate  type 
as  opposed  to  solid  type  elements. 

On  the  basis  of  these  requirements  and  considering  viable 
alternatives,  three  suitable  continuum  formulations  have  been 
developed  and  are  herein  denoted  as  the  (i)  Higher  Order  Displacement, 
(ii)  Modi fied-Kirchhoff  and  (iii)  Hybrid  Stress  formulations, 
respectively.  The  former  two  formulations  have  been  incorporated  in  a 
computer  code  and  the  various  elements  have  been  tested  on  the  basis 
of  correlations  with  known  analytical,  numerical,  and  experimental 
solutions.  Numerous  tests  have  been  performed  for  linear  static  and 
linear  dynamic  cases.  It  is  noted  that  the  code  has  some  unique 
features,  e.g.,  it  can  assemble  elements  having  an  unequal  number  of 
degrees  of  freedom  at  its  nodes,  it  treats  arbitrary  ply  orientations 
and  it  performs  integration  on  a  layer-by-layer  basis  through  the 
laminate.  Herein  a  layer  refers  to  either  a  lamina  or  to  a  sub-set  of 
laminae  having  equal  ply  orientations.  The  latter  feature  is 
essential  in  developing  a  fully  nonlinear  capability. 

Significant  efforts  have  also  been  devoted  to  developing  a 
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suitable  large  displacement  formulation.  Due  to  the  requirement  that 
interlaminar  stresses  be  accurately  represented,  a  total  Lagrangian 
formulation  is  utilized  and  is  based  upon  the  complete  Green's  strain 
tensor.  A  geometric  and  large-displacement  stiffness  formulation  has 
been  implemented  in  the  computer  code  based  upon  a  form  of  the 
nonlinear  strain-nodal  displacement  relationships  suitable  for  each  of 
the  elements  under  development. 

An  extensive  literature  survey  has  been  performed  to  identify 
analytically  tractable  methods  of  treating  damage  accumulation  in 
composites.  Since  emphasis  in  this  work  is  on  the  development  of 
incrementil  response  solutions,  the  computational  approach  must  have 
the  capability  to  (i)  predict  and  differentiate  between  relevant 
failure  inodes,  (ii)  modify  constitutive  equations  appropriately  and 
(iii)  perform  equilibrium  iterations  to  assure  stress  redistribution 
based  upon  the  extent  of  damage.  Use  of  "piecewise  smooth"  failure 
cri  teria  in  conjunction  with  "damage  state"  variables  provide  a  good 
basis  for  i ncremental ly  tracking  damage.  This  approach  has  been 
incorporated  in  the  computer  code.  Note  that  integration  for  an 
element  is  performed  on  a  layer-by-layer  basis  which  allows  for  damage 
effects  to  be  characterized  at  the  layer  level. 

Experimental  data  of  the  type  required  to  substantiate  danvage 
predictions  has  been  assembled  to  the  extent  possible.  Analysis/test 
correlations  have  been  performed  for  selected  laminates.  It  is  noted 
that  useful  experimental  data  is  quite  limited. 

Technical  progress  in  this  program  has  been  substantially  on 
schedule.  It  is  believed  that  the  originally  proposed  three  year 
program  can  be  completed  within  the  given  time  frame. 


II.  SUMMARY  BY  TASK 


This  section  presents  technical  highlights  of  the  research 
efforts  to  date  for  each  of  the  three  tasks.  Details  of  the 
analytical  formulation  are  presented  in  the  Appendices. 


II.l.  TASK  I:  Nonlinear  Displacement  Formulation  for  Composite  Media 


II. 1.1  Continuum  Formulation 


Two  variational  principles,  the  principle  of  Minimum  Potential 
Energy  and  the  Principle  of  Modified  Complementary  Energy,  are  used  to 
develop  t</o  distinctly  different  finite  element  models,  the  assumed 
displacement  model  and  the  hybrid  stress  model  respectively.  These 
models  incorporate  the  effects  of  transverse  shear  and  normal 
deformations  whose  contributions  are  recognized  as  essential  for 
accurate  laminate  analysis  [1-10]. 

Within  each  formulation,  element  stiffness  and  force  matrices  are 
determined  for  each  element,  these  matrices  are  then  assembled  to 
represent  the  final  system  of  equations  and  a  solution  procedure  for 
tile  unknown  nodal  displacements  is  provided.  Coordinate 
transformations  to  describe  ply  orientations  of  a  composite  media  are 
taken  into  account.  The  in-plane  stresses  are  calculated  from 
constitutive  relations  of  ortiiotropic  continuum  whereas  transverse 
shear  and  normal  stresses  are  calculated  from  equilibrium 
considerations.  At  present,  emphasis  is  placed  on  the  displacement 
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based  iTTOdols  and  these  have  been  tested  for  linear  static  and  dynamic 
analysis.  The  test  problems  and  the  results  are  presented  in  Section 
II. 1.4.  The  finite  element  models  are  herein  briefly  discussed. 

i.  Assumed  Displacement  Model 

A.  Higher  Order  Displacement  Formulation 

The  thru-the-thickness  effects  can  be  incorporated  into  an 
analysis  by  choosing  a  displacement  field  that  eliminates  two  n>ajor 
shortcomings  of  the  classical  plate  theory;  namely  normals  remain 
normal  and  in-plane  displacements  are  linear  thru  the  thickness. 

These  shortcomings  are  eliminated  by  prescribing  independently  the 
reference  surface  displacements  and  rotations  of  the  normal  and 
including  higher  order  terms  for  in-plane  displacements.  This  is 
a ccompl  i '.hod  by  the  following  variation 

-.y,-')  =  UQ(;<,y)  +  z.j((x,y)  +  z'  •x(x,y) 

'.-w'.  ,y,z)  =  Vo(x,y)  +  z.,,(x,y)  +  z-:y{x,y) 

■'d  /  ,  y ,  z  )  =  v/g  ( X  ,  y ) 


The  neutral  surface  di solacements  are  represented  by  u  ,  v  and  w  , 

o  c 

the  rotation  about  y-axis  is  denoted  by  and  the  rotation  about  the 
x-axis  is  x,,.  The  coefficients  of  z  ,  i.e.,  and  s‘y»  ^*^2 
contributions  from  transverse  deformations  [5,6]. 

The  elements  developed  are  designated  as  the  quadrilateral  higher 
order  displacement  (QHD)  models.  QHD40  is  an  eight-noded  element  with 
seven  degrees  of  freedom  (three  midsurface  displacements,  two 


rotations  and  tvo  higher  order  terms  tor  in-plane  displacements)  per 
corner  node  and  three  degrees  of  freedom  {transverse  midsurface 
displacement  and  two  rotations)  per  mid-side  node.  Element  QHD28  is 
a  simplified  version  of  QHD40  where  the  mid-side  nodes  are  eliminated. 

It  should  be  noted  that  when  the  tivo  higher  order  terms  for  in-plane 
displacements  at  each  corner  node  are  omitted,  QHD28  reduces  to  the 
widely  used  four-noded  bilinear  plate  element  (0HD2O). 

The  transverse  shear  and  normal  stresses  of  QHD40  display  a  cubic 
variation  thru-the-thickness.  The  displacement  field,  nodal  degrees 
of  freedom  and  tiie  resulting  stress  fields  are  stated  in  Appendix  lA. 

B.  Modi fied-Ki rchhof f  Formulation 

The  Ki ''chhof f-Love  assumption  for  normals  to  the  reference 
surface  is  relaxed  by  incorporating  shear  rotations  as  additional 
degrees  of  freedom  in  the  formulation  [10].  Thus  the  assumed 
displacement  field  allows  the  transverse  shear  Jeforina tions  but 
neglijcts  the  transverse  normal  deformia  ti  ans.  The  rotations  T,.and  are 
incorporated  in  the  displacement  variation  as  follows 

w(x,y)  =  Wo(x,y) 
u(x,y,z)  =  iio(x,y)  -  z{  ~+ 

X 

v(x,y.z)  -  Vo(x,y)  -  ‘  '-y) 

The  trmsverse  displacement  w(x,y)  is  chosen  such  tiiat  it  will 
guarantee  plausible  stress  fields  which  will  characterize  the 
trinsversi?  effects  accurately. 


This  approach  is  implemented  in  the  formulation  of  an  eight-node 
quadrilateral  element  with  32  degrees  of  freedom-  QD32,  a  six-node 
triangular  element  with  27  d.o.f.-  TD27  and  a  seven-node  triangulat 
element  with  27  d.o.f.-  TD27M.  The  stress  fields  obtained  for  these 
elements  represents  a  quadratic  thru  the  thickness  variation  for  the 
transverse  shear  stresses  and  a  cubic  variation  for  the  transverse 
normal  stress.  The  respective  displacement  fields,  nodal  degrees  of 
freedom  and  stress  fields  are  given  in  Appendix  IB. 

ii.  Hybrid  Stress  Model 

In  this  formulation  a  stress  distribution  within  the  interior  of 
the  element  is  expressed  in  terms  of  finite  parameters  such  that 
equilibrium  is  satisfied,  also  an  assumed  displacement  distribution  is 
used  on  the  boundary  of  the  element  expressed  in  terms  of  generalized 
nodal  displacements  such  that  the  interelement  compatibility  is 
retained  [4]. 

The  element  developed,  QHS32,  is  a  four-node  quadrilateral  with 
32  .legrees  of  freedom.  In  addition  to  an  assumed  di  splacement  field 
it  dis  a  2o-parame ter  stress  field  which  provides  cubic  variation  for 
transverse  siiear  stresses  and  a  quartic  variation  for  transverse 
nor'i.ial  stress  timough  tlie  thickness  of  the  laminate.  The  stress  field 
O  3ny  wi  t:i  th'’  issumed  displacement  variation  is  st.itod  in  Appendix  II 

1I-L.2.  l.arge  dripiucement  i-ormulation 

Inclusion  of  jeometri cal ly  nonlinear  effects  in  the  formulation 
must  be  based  upon  both  the  geomietry  to  be  analyzed  and  upon  the  type 
of  stress  prediction  capabilities  desired.  The  classical  approach  to 
thin  plate  analysis  has  been  to  use  the  Kirchhof f-Love  assumptions  in 


conjunction  with  the  nonlinear  von  Karman  relations  [11,12],  As 
previously  indicated,  the  Ki rchhof f-Love  assumptions  are  relaxed  in 
this  work  to  allow  for  a  more  accurate  definition  of 
interlaminar-shear  and  interlaminar-normal  stress  variations.  Thes 
stresses  can  vary  substantially  through-the-thickness  for  the 
geometries  of  interest,  i.e.,  thin  to  moderately  thick  plate  type 
structures.  Furthermore,  the  requirement  that  these  stresses  be 
accurately  determined  means  that  the  nonlinear  portion  of  the 
strain-di splacement  relationship  must  contain  all  significant 
coordinate  displacements.  The  complete  Green's  strain  tensor  is 
utilized  in  this  work,  therefore,  to  account  for  all  significant 
contributions  to  the  interlaminar  stress  field.  Wi tn  respect  to 
fixed  Cartesian  coordinates,  x,  y,  and  z,  the  strain  tensor  has  the 
f  LU'in 


X 


■■  / 


■  w 


where  u,  v  and  w  represent  displacements  inb  the  x,y,z  coordinate 
directions,  respectively.  Note  that  the  other  strain  components  ar 
oijtained  by  a  suitable  permuta  tion.  In  srnal  1-di  splacernent  analysis 
tne  (juidratic  terms  are  neglected  to  give  simply  the  linear  strain 
i  ppi'ox  i  mi  Li  on . 

liised  on  the  Green's  strain  tensor,  the  strain  to  nodal  point 


1  i  spl  i  ^  ■  mol  t  relationshif'  can  be  specified  tor  elements  inder 


J -'volopinent.  It  tiki^s  th-:?  form 


wher'j  II  is  the  vector  of  strain  components,  |^'.|  the  vector  of 
nodal  point  displacements  and  [3]  a  function  of  derivatives  of  the 
element  shape  functions.  The  quadratic  terms  in  the  strain  tensor 
result  in  [B]  being  a  function  of  displacement  state  and,  therefore, 
an  incremental  equilibrium  formulation  is  required.  The  incremental 
strain-nodal  di  sj)!  acement  relationship  takes  the  form 

=  (i3p  .  IBJI 

where  {  |  and  |  |  represent  incremental  strains  and  nodal 

displacements,  respectively,  [3  ]  and  [Bj]  are  the  small  and  large 
displacement  contributions  to  the  incremental  strains.  Based  on  the 
incremental  i;gui  1  i orium  equations,  the  displacement  formulation  gives 
the  i'orce-di  .placenont  r ila ti onshi ps 


"L 


where  [0]  is  in  elisticity  matrix  obtained  simply  from  the 
c onsti  tuti vj  iquitiins  and  integration  is  over  tile  volume  V  of  the 
element.  [K  J  is  lenoihni  tiie  sma  1 1 -di splacern'cn t  stiffness  'matrix  and 
[<jl  is  I'Mioted  th:  1  .i  rgo-d  i  spl 'ic  amen  t  stiffness  matrix.  Since 
■^esponse  is  ilso  .i  function  of  stress  stite,  the  geometrical  stiffness 
es  rv'jU  L  ri-'d  intl  is 


i  r.r  I  < 


I'bi  iinod  from 


where  |a|  is  the  vector  of  stress  components.  Note  that  the  hybrid 
stress  formulation  similarly  gives  force-displacement  forms  involving 
the  stress  and  displacement  state. 

Inertial  effects  are  analytically  treated  as  a  mass  matrix  Cm] 
which  is  a  function  of  density  and  the  element  shape  functions  (see 
Appendix  III).  These  matrix  forms  are  required  in  formulating 
static/dynamic  response  solutions  and  the  incremental  equilibrium  have 
the  general  form 

(.MKcu)  +  (iKj  t  fKj  .  ^ 

where  the  rass  and  stiffness  matrices  represent  an  assembly  of  the 
elemental  matrices  previously  discussed,  |-'u|  and  |<5u[  represent  the 
incremental  displacements  and  accelerations  for  the  mathematical  model 
and  I'l'}  represents  the  vector  of  incrementally  applied  forces. 

In  developing  a  geometrically  nonlinear  formulation,  the  effort 
is  largely  in  defining  the  incremental  strain-nodal  displacement 
relationship.  Having  developed  this  relationship  for  a  particular 
element,  stiffness  irva trices  are  readily  developed  as  the  preceding 
equations  indicati?.  These  relationships  are  presented  in  Appendix  IV. 

The  form  of  these  equations  is  the  same  for  all  elements. 

II. 1.3.  Computer  Implementation 

A  computer  code  has  been  developed  for  the  purpose  of 
impleineiiting  the  various  continuum  formulations.  At  present,  the  code 
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performs  the  following  functions: 

element  stiffness  matrix  (linear,  nonlinear,  geometric) 
generation 

element  mass  matrix  generation 
assembly  of  equilibrium  equations 
decomposition  and  solution  of  equilibrium  equations 
equilibrium  iteration  for  incremental  solutions 
fundamental  frequency  and  mode  shape  calculation 
A  characteristic  of  the  elements  under  development  is  that  node 
points  can  have  different  numbers  of  degrees  of  freedom,  i.e., 
typically  mid-side  nodes  have  fewer  degrees  of  freedom  than  corner 
nodes.  The  code  has  been  fashioned  to  handle  this  condition.  All  of 
the  integration  is  performed  on  a  layer-by-layer  basis  thru  the 
thickness  of  the  laminate.  This  approach  is  fundamental  to  developing 
the  capability  to  allow  for  inelastic  material  behavior  and, 
ultimately,  to  tile  inclusion  of  damage  mechanisms  in  the  formulation. 

Since  solution  of  the  equilibrium  equations  is  a  vital  component 
in  the  overall  solution  strategy,  it  is  appropriate  to  discuss  the 
numerical  methodology  used  in  solving  these  equations.  The  intent  is 
to  obtain  a  higher  ordered  variation  of  the  transverse  shear  and 
normal  stresses  (  *  ,  a  ,  and  a  )  tnan  can  be  obtained  via  the 

X  z  y  z  z  z 

c.institut. ive  equations.  The  solution  procedure  can  be  thought  of  as 
described  below.  Assume  that  the  in-plane  stresses  ) 

within  each  layer  of  a  particular  element  have  been  determined  at 
selected  locations,  i.e.,  through  solution  of  the  constitutive 
equations.  In  the  code  as  presently  written,  these  locations  are 
specified  as  the  element  centroid  and  element  nodal  points.  The 
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equilibrium  equations  (in  the  absence  of  body  forces)  have  the 
indicia!  form 


^ij,j  =  0 


from  which  it  follows  that  the  thru-the-thickness  shear  stress 
variation  can  be  written  in  numerical  form  for  the  i*^^  layer  as 


"xzi  -  -(axx.x  +  -z,. 


"  '^'xy,x  +  -yy,y)^  .I2i 


Here,  the  left-hand-side  represents  the  change  in  stress  from  the 
lower  to  the  upper  surface  of  the  i*^^layer  and  is  the  thickness  of 
the  i’^*'layer  at  a  particular  location.  The  derivatives  with  respect 
to  X  and  y  in  the  expressions  above  are  readily  computed;  this  is 
because  in-plane  stresses  within  a  layer  are  related  to  element 
displacements  through  derivatives  of  element  shape  functions  in 
conjunction  with  a  riiatGridl  definition. 

For  an  n  layered  laminate,  n  equations  can  be  written  in  terms  of 
both  the  unknown  shear  stresses  at  layer  interfaces  and  the  shear 
stresses  at  tlie  laminate  surfaces.  Assuming  the  laminate  has 
shear-free  surfaces,  the  equations  above  give  n  equations  in  n-1 
unknowns,  so  that,  the  equation  set  is  over-determined.  The  equations 
have  the  matrix  form  below 


C  V'  v'-T 


where  I^zi"  ~  ^'^xx,x  ^  '^xy,y)^  and  ^^®Pi^®sents  the  shear  stress 

acting  at  the  interface  of  the  j-1 and  layer.  A  similar 
equation  set  is  obtained  by  replacing  cr^^j  ^yzj  ^xzi  ^yzi, 

These  equations  are  solved  by  utilizing  a  least-squares 
orthonormalization  procedure  [13],  Due  to  the  simplicity  of  the  terms 
in  the  coefficient  matrix,  a  concise  closed-form  solution  is  obtained. 
Having  determined  the  transverse  shear  stresses,  the  transverse 
normal  stress  variation  is  determined  through  the  numerical  form  of 
the  third  equilibrium  equation  for  the  i  layer 


'xz’x  ^’y2»y^-j 


As  before,  the  left-hand-side  represents  the  change  in  stress  through 

til 

tile  i  layer.  Appropriate  polynomial  functions  are  utilized  to 
describe  the  c  .  and  a.,„  in-plane  variation.  These  functions  are 
differentiated  to  obtain  the  right-hand-side  of  the  equations  above. 
Again  the  equation  set  is  overdetermined  because  the  normal  tractions 
are  known  at  the  laminate  surfaces.  Solving  fora22  proceeds, 
therefore,  in  identically  the  same  nvanner  as  discussed  in  calculating 
and  Parenthetically,  inclusion  of  body  forces  at  a  later  date 


1 

1  1 
-1  1 


•u  •**,*  ■ 

.  • .  • . 


^  -'v'' 

•  •  •  -  o  •  • 

.ywV.'. 


n  X  (n-1)  (n-1)  x  1 


(n  X  1 ) 


#  . 


can  be  accomplished  with  little  difficulty. 


It  should  be  emphasized  that,  the  successful  application  of 
Higher  Order  Displacement  type  elements,  i.e.,  for  particularly  thin 
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geometries,  i s  to  utilize  reduced  numerical  integration  where  as  this 
is  not  necessary  for  the  Modified  Kirchhoff  formulation.  This 
approximation  technique  brings  along  the  choice  of  implementing  it 
overall  or  selectively  to  the  strain  energy  components.  For  the  QHD 
formulation,  only  the  transverse  shear  components  are  integrated  with 
reduced  order  [14-16].  An  undesirable  aspect  of  this  approach  is  that 
the  reduced  integration  order  may  affect  the  physical  behavior  of  the 
element  by  introducing  spurious  zero  energy  modes.  It  is  desirable  to 
have  only  rigid  body  modes  since  there  does  not  yet  seem  to  be  a 
generally  accepted  method  of  controlling  the  additional  modes. 


1 1.1.4.  Analytical  Verification 

As  noted,  elements  formulated  on  the  basis  of  independent 
transverse  di splacements  and  rotations,  require  reduced  quadrature  for 
good  performance.  For  QHD40,  3x3  Gaussian  quadrature  along  with  the 
2x2  quadrature  for  tlie  transverse  shear  components  in  employed.  QH028 
and  QHL)20  formulations  are  similarly  integrated  with  2x2  and  1x1 
Gaussian  quadratures.  Manipulation  of  quadrature  rules  may  produce 
spurious  zero  energy  modes  in  addition  to  the  required  rigid  body 
modes,  thus  detracting  from  overall  element  performance.  [16,17].  A 
spectral  (eigenvalue)  test  has  been  conducted  with  and  without  full 
quadrature  to  observe  the  zero  energy  modes  of  the  QHD  elements.  The 
quadrature  order,  the  number  of  zero  eigenvalues  and  the  corresponding 
number  of  spurious  zero  energy  modes  for  the  QHD40,  QHD28,  and  QH020 
elements  are  listed  in  Table  1.  The  spurious  mode  shapes  associated 


with  tile  QH028  element  are  illustrated  in  Figure  1.  Since  the  QHD40 
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formulation  does  not  exhibit  spurious  modes,  it  can  be  utilized  in 
modelling  complex  geometries  without  concern  for  controlling  such 
behavi or. 

It  is  also  noteworthy  to  observe  the  effect  of  reduced 
integration  on  the  representation  of  the  generalized  forces.  In 
order  to  illustrate  the  effect,  the  forces  associated  with  the 
transverse  displacement  of  a  corner  node  are  sketched  in  Fig.  2  for 
QHD28  with  and  without  reduced  integration  respectively. 

In  the  examples  that  follow,  performance  of  the  QHD  formulation 
is  demonstrated  by  comparing  results  to  those  obtained  by  classical 
plate  theory  (CRT),  by  elasticity  and  by  other  finite  element 
formulations  for  linear  static  and  linear  dynamic  analyses.  Some 
results  are  also  presented  for  the  QD  formulation.  The  orthotropic 
material  properties  used  throughout  are  tabulated  in  Table  2. 
Geometries  studied  include  cylindrical  bending  of  a  plate  as  well  as 
bending  of  simply  supported  square  and  rectangular  plates.  Various 
ply  layups  are  considered  and  loading  is  that  of  a  sinusoidally  and 
uniformly  distributed  pressure.  Cylindrical  bending  is  modelled  as  a 
strip  of  twenty  elements.  For  the  simply  supported  plates,  symmetry 
considerations  allow  that  only  a  quadrant  of  the  plate  need  be 
modelled.  Fineness  of  the  mesh  is  varied  to  demonstrate  solution 
convergence.  Additionally,  distorted  meshes  are  considered  to 
demonstrate  inodelling  considerations.  For  the  examples  involving 
symmetric  layups,  the  quadratic  terms  of  QHU40  are  restrained;  so 
that,  32  degree  of  freedom  elements  are  utilized  to  ohtciin  these 
solutions.  This  is  allowable  in  these  particular  cases  because  the 
quadratic  terms  do  not  significantly  affect  the  results.  This  is  not 


true  in  the  first  example  considered. 

Statio  Response  CalauZationc 

Cylindrical  Bending-  Bidirectional  (0  /90)  Sine  Load,  Material  II 

Fibers  run  parallel  to  the  plane  of  curvature  in  the  lower  layer 
and  are  rotated  90°  in  the  upper  layer  of  the  plate.  Layers  are  of 
equal  thickness  which  is  true  in  the  subsequent  example  problems  as 
well.  The  elasticity  solution  obtained  by  Pagano  [9]  gives  a  nearly 
quadratic  z  variation  in  u,  where  u  is  the  normalized  in-plane 
displacement  of  the  laminate.  In  this  instance,  inclusion  of  the  z“ 
terms  in  the  finite  element  modelling  should  affect  the  results.  This 
is  demonstrated  in  Figures  3  to  5.  Results  demonstrate  differences 
obtai.iod  with  and  without  quadratic  terms.  The  difference  is  greatest 
for  the  lower  aspect  ratios,  e.g.,  for  S  =  4  a  difference  of  12%  is 
obtained.  In  Figure  4,  the  calculated  normalized  in-plane  stress 
variation  is  presented  for  an  aspect  ratio  of  4.  Note  that  maximum 
variation  is  presented  for  an  aspect  ratio  of  4.  Note  that  maximum 
stresses  differ  by  some  36?  when  computed  with  and  without  the  z 
terms,  respectively.  The  effect  of  including  quadratic  terms  in  the 
finite  element  solution  is,  therefore,  much  more  pronounced  when 
stresses  as  opposed  to  displacements  are  considered.  Figure  5 
demonstrates  this  effect  on  stress  computation  as  a  function  of  aspect 
ra.io.  Note  that  calculated  quantities  are  normalized  in  this  example 
and  in  those  that  follow  as  in  the  cited  references. 

Cylindrical  Bending-Symmetric  (0  /90  /0)Sine  Load,  Material  II 

For  this  geometry,  fibers  are  parallel  to  the  plane  of  curvature 


in  the  outer  layers  and  rotated  90°  in  the  middle  layer.  Calculated 
stresses  are  compared  to  the  elasticity  solution  of  Pagano  [9]. 

Figures  6  and  7  present  the  normalized  transverse  shear  stress 
variation  at  the  simply  supported  boundary.  Figures  3  and  9  present 
the  normalized  in-plane  stress  variation  at  the  center  of  the  bent 
surface. 

Simply  Supported  Square  Plate-  (0  /90  /O)  Sine  Load,  Material  II 

Fibers  in  the  outer  layers  of  the  laminate  run  parallel  to  the  x 
axis  while  those  in  the  middle  layer  run  parallel  to  the  y  axis,  where 
the  origin  of  coordinates  is  located  at  a  corner  of  the  plate  and  in 
the  mid-plane  (see  Figure  10).  This  coordinate  system  is  consistent 
with  examples  that  follow  as  well.  Consider  the  plate  as  having 
planar  dimensions  a  x  a  and  total  thickness  h.  Solutions  have  been 
generated  for  aspect  ratios  S  =  4  to  100,  where  S  =  a/h.  Transverse 
shear  stress  variation  at  (x,y)  coordinates  {0,a/2)  and  in-plane 
shear  stress  variation  at  coordinates  (0,0)  are  presented  in 
Figures  10  and  11  for  an  aspect  ratio  of  4.  Note  that  the  comparison 
is  between  the  present  finite  element  results  and  those  obtained  via 
elasticity  [18]  and  CPT.  Calculated  short-transverse  normal  stress 
variation  is  presented  in  Figures  12  and  13.  These  stresses  are 
normalized  as  0^2=  at  the  center  of  the  plate  and  as 

at  the  edge  of  the  plate.  Results  are  compared  to  those  obtained  by 
elasticity  over  a  range  of  aspect  ratios  in  Table  3.  Similar  results 
are  given  in  Table  3.1  for  the  QO  formulation.  Convergence 
characteristics  are  demonstrated  by  presenting  results  obtained  using 
2x2,  3x3,  and  6x6  meshes.  The  finer  mesh  gives  better  agreement,  but 
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the  coarser  mesh  gives  very  reasonable  correlation  also. 

The  effects  of  distorting  the  mesh  have  also  been  considered  to  a 
limited  extent.  Results  have  been  obtained  for  the  relatively  coarse 
meshes  shown  in  Figure  14.  Calculated  stresses  and  displacements  are 
presented  as  a  function  of  aspect  ratio  and  compared  to  the  elasticity 
solutins  in  Table  4.  As  expected,  the  values  are  not  as  accurately 
determined  as  are  those  obtained  via  the  regular  meshes.  Distortion 
of  the  mesh  has  a  much  more  dramatic  effect  upon  the  calculated 
transverse  shear  stresses  than  upon  the  calculated  in-plane  stresses 
and  displacements.  Since  the  transverse  stresses  are  based  on 
equilibrium  considerations,  it  seems  the  mesh  must  be  refined  enough 
to  reasonably  approximate  equilibrium.  This  is  especially  apparent  in 
comparing  results  obtained  for  mesh  A  to  those  obtained  for  mesh  C. 

In  each  of  these  cases,  elements  having  a  taper  ratio  of  2  to  1  are 
utilized.  Mesh  C  gives  significantly  improved  transverse  stresses, 
hov^ever,  because  the  mesh  is  fine  enough  to  better  represent  the 
loading  distribution. 


Simply  Supported  Rectangular  Plate(0  /90  /O)  Sine  Load,  Material  II 
Orthotropic  layers  have  the  same  orientation  as  in  the  previous 
example.  The  plate  has  dimensions  a  x  b,  where  b  is  three  times  a. 
Solutions  have  been  obtained  tor  aspect  ratios  (S  =  a/h)  ranging  from 
4  to  100.  Transverse  shear  stress  variation  at  coordinates  (a/2,0) 
is  given  in  Figure  15  tor  an  aspect  ratio  S  =  4.  Comparison  is  made 
to  both  elasticity  and  CPT  solutions.  A  full  range  of  results  are 
presented  in  Table  5  and  compared  to  those  obtained  via  elasticity 
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[18]  and  to  those  obtained  by  Reddy  [19]  in  a  recent  finite  element 
formulation.  Correlation  with  elasticity  is  quite  good,  particularly 
for  aspect  ratios  of  10  and  above,  and  appear  to  be  nrora  accurate  than 
those  obtained  with  the  alternate  finite  element  solution. 

Simply  Supported  Square  Plate  (0  /90  /90  /O)  Sine  Load,  Material  II 

The  laminate  geometry  consists  of  outer  layers  with  fibers 
parallel  to  the  x  axis  and  inner  layers  with  fibers  parallel  to  the  y 
axis.  The  plate  has  planar  dimension  a  x  a  and  total  thickness  h. 
Stress  and  displacement  results  are  presented  in  Table  6  for  aspect 
ratios  ranging  from  4  to  100.  Similar  results  are  given  in  Table  6.1 
based  on  the  QU  formulation.  Results  are  compared  to  both  elasticity 
and  to  other  finite  element  results.  Again,  the  computed  values  are 
in  excellent  agreement  with  elasticity  [20]  for  moderately  thick  to 
thin  geometries  and  are  more  accurate  than,  the  compared  to  numerical 
resul ts. 

Solutions  also  have  been  obtained  for  the  present  geometry  on  the 
basis  of  reduced  vs.  full  integration.  This  comparison  is 
demonstrated  in  Figures  lo  and  17  by  giving  percent  error  in 
calculated  values  vs.  aspect  ratio.  It  is  apparent  that  reduced 
integration  is  particularly  needed  to  minimize  errors  in  calculated 
transverse  stresses  and,  furthermore,  solution  validity  over  a  wide 
range  of  laminate  geometries  is  demonstrated. 

To  assess  the  effects  of  finite  element  formulations,  aspect 
ratio,  support  conditions  and  the  lamina  stacking  sequences  on  the 
fundamental  natural  frequencies  of  composite  plates,  the  problems 


listsa  in  Taoie  7  are  considered  [21] 

The  non-dimensional ized  fundamental  frequency  for  the  cross-ply 
laminate  of  ProDlem  1  versus  aspect  ratios  is  given  in  Table  8.  As 
can  be  seen,  all  three  elements  predict  frequencies  that  are  in 
excellent  agreement  with  the  closed  form  solutions  obtained  by  Reddy 
[22J. 

The  effects  of  higher  order  terms  in  the  displacement  based 
finite  element  formulations  are  investigated  for  ProDlem  2.  Here,  the 
performances  of  QHD40  and  QHD23  (with  higher  order  terms  locked)  are 
compared  to  elements  STPOl  and  STP03  of  [23]  with  linear  and  cubic 
variations  through  tiie  thickness  respectively.  The  results  are 
summarized  in  Table  9.  The  normalized  fundamental  frequencies  of 
Problem  3  are  displayed  in  Figure  13.  Note  that  the 
non-dimensional ized  fundamentil  frequency  increases  as  the  angle  of 
orientation  is  increased  for  both  symmetric  and  antisymmetric 
ang1e-ply  square  plates.  This  observation  is  in  excellent  agreement 
with  Reddy's  [22]  antisymmetric  laminate.  In  Figure  19,  a  decrease  in 
the  fundamental  f'-equency  is  observed  as  the  angle  of  orientation  is 
increased  for  the  angle-ply,  cantilever,  rectangular  and  square  plates 
of  Problem  4.  The  difference  between  Figures  13  and  19  are  attributed 
primarily  to  the  different  support  conditions. 

Further  investigations,  Problein  5,  of  angle-ply  laminates  are 
summari zed  in  Table  10.  The  stacking  sequences  of  reference  [24]  are 
used  to  illustrate  their  effects  on  the  fundamental  frequency 
calcula tions.  The  numbers  within  parenthesis  are  calculated  by 
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Shown  in  Figure  20  is  the  variation  of  the  non-dimensional ized 
fundamental  frequency  for  cylindrical  bending  problem,  as  calculated 
via  the  QD  formulation  and  the  classical  plate  theory.  For  comparison 
purposes,  the  frequencies  are  normalized  with  respect  to  the  classical 
plate  theory  results. 

Ti'aKcicKt  Ft  j?  cr.r-j  Calculationc 

Element  performance  has  been  evaluated  with  respect  to  predicting 
linear-transient  response.  Both  displacements  and  stresses  have  been 
determined  for  a  variety  of  laminated  plate  geometries  subjected  to 
instantaneously  applied  pressure  loading.  These  results  have  been 
compared  to  those  obtained  via  both  CRT  and  a  shear  deformable  theory 
(SOT)  [26].  Typical  results  are  presented  in  Figures  21  and  22  for  a 
(0/90)  square  plate  [27].  In  this  example,  the  plate  is  quite  thick 
in  that  it  has  an  aspect  ratio  of  5. 

II. 2.  TASK  II;  Incorporate  Damage  Mechanisms  into  Dynamic  Response 
Formulation 

The  literature  survey  [23-63]  perfonned  has  been  quite  helpful  in 
ter:ns  of  delineating  the  viable  approaches  to  including  damage 
meclianisms  in  the  analysis.  Relevant  failure  modes  of  interest 
include  tiios*;  listed  below 

(i)  fiber  fracture 

(ii)  fiuer-inatrix  debonding 

(iii)  matrix  cracking  (parallel  and  transverse  to  fibers) 

(iv)  delamination 

Several  snooth  fiilure  criteria,  e.g.,  [64-67]  have  been  developed  in 
recent  years  to  represent  the  failure  of  composites.  These  criteria, 
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to  varying  ueyrees,  can  predict  "failure"  but  do  not  identify  a 
particular  mode  of  failure.  In  performing  incremental  "damage" 
analysis,  it  is  essential  to  both  predict  failure  and  to  characterize 
it,  e.g.,  do  fibers  rupture,  does  delamination  occur,  etc.  The 
computational  approach  must,  therefore,  differentiate  between  viable 
failure  modes  and  appropriately  alter  the  constitutive  equations  on 
an  incremental  basis.  This  can  be  accomplished  by  impementing  a 
piecewise  smooth  failure  criteria,  e.g.,  [28]  in  the  finite  element 
formulation.  The  general  failure  criteria  is  then  comprised  of  m 
separate  inequalities  of  the  form 

at  the  layer  level  within  each  element.  These  criteria  can 
differentiate  between  (i)  tensile  and  compressive  fiber  failure,  (ii) 
tensile  ind  compressive  matrix  failure  and  (iii)  delamination  at 
layer  interfaces. 

As  progressive  damage  occurs  throughout  incremental  loading 
(whether)  it  be  static  or  dynamic),  it  is  essential  that  violation  of 
failure  criteria  inequalities  be  reflected  in  modification  of  the 
material  properties.  This  can  be  achieved  by  including  damage  state 
variables  [47]  in  the  constitutive  equations  to  reflect  "stiffness 
reduction."  These  equations  can  be  represented  as 

■  ^  rOjL:];-  ■ 

where  [D]  represijnts  tiie  laiterial  matrix  and  [  ■  ]  contains  the  damage 
state  variables.  The  latter  provide  the  basis  for  changing  the  Ojj 
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terms  based  upon  the  extent  to  which  the  failure  criteria  are 
vi olated. 

In  conjunction  with  the  above  it  is  essential  to  perform 
equilibrium  iterations  within  each  analysis  increment.  This  is 
required  to  assure  that  stress  redistribution  is  properly  accounted 
for  as  damage  progresses. 

The  relevant  failure  nodes  of  interest  and  appropriate  criteria 
used  in  the  incremental  analysis  are  listed  in  Table  II. 


The  damage  histories  of  graphite/epoxy  {T30b/5208)  composite 
Icami  nates  under  in-plane  load  increments  are  presented  for  selected 
models  [68],  The  ability  to  differentiate  between  relevant 
failure/damage  modes  is  illustrated  in  Figures  23  and  24. 

11,3.3  TASK  III:  Correlation  of  Formulated  Response  Model  with 
Experimental  Data 

Some  quantitative  data  relating  to  the  impact  damage  of  composite 
specimens  has  been  assembled  [69-76].  It  will  be  utilized  along  with 
any  additional  daU  obtained  to  perform  analysi s/tes t  correlations. 
Since  the  nonlinear  formulation  including  damage  effects  is  not 
complete,  no  use  of  the  test  data  has  been  made  to  this  point. 
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Table  1.  Spurious  Zero  Energy  Modes  of  the  QHD  Family 


Quadrature  Order 


Number  of  Zero 
Eigenvalues 


Number  of  Spurious 


QHD40 

— 

3x3  with  2x2 
for  transverse 
shear  terms 

6 

0 

QHD28 

2x2  with  1x1 
for  transverse 
shear  terms 

9 

3 

QHD20 

2x2  with  1x1 
for  transverse 
shear  terms 

8 

2 

34 


CPI  .534  .IKO  .0213  .395  .0823 


35 


TABLE  3.1 


Approach 


FEM  (3x3  Mesh) 
FEM  (6x6  Mesh) 
Elasticity 


X 

fil.  +!!■) 
'2* 


.399 

.392 

.755 


^2’  2’  ~6 


.562 

.543 

.556 


XV 


(0.  0,  ±-^') 


.0513 

.0463 

.0505 


(0,  0) 


.372 

.357 

.282 


y  z 


(^.  0,  0) 


.304 

.280 

.217 


10 


FEM  (3x3  Mesh) 
FEM  (6x6  Mesh) 
Elasticity 


.514 

.502 

.590 


.246 

.270 

.288 


.0299 

.0284 

.0289 


.406 

.387 

.357 


.175 

.142 

.123 


20 


50 


FEM  (3x3  Mesh) 
FEM  (6x6  Mesh 
Elasticity 


FEM  (3x3  Mesh) 
FEM  (6x6  Mesh) 
Elasticity 


.547 

,533 

.552 


.558 

,543 

.541 


.157 

.186 

.210 


.0245 

.0234 

.0234 


.128 

.159 

.185 


.0227 

.0219 

.0216 


.418 

.398 

.385 


.423 

.402 

.393 


.141 

.107 

.0938 


.130 

.0961 

.0842 


100 


FEM  (3x3  Mesh) 
FEM  (6x6  Mesh) 
Elasticity 


.359 

.544 

,539 


.123 

.155 

.181 


.0225 

.0216 

.0213 


.423 

.403 

.395 


.128 

.0944 

.0828 


■/O'J 


•5  39  .IKl  .02  It  .2<r) 


TABLE  5 

STRESSES  AND  DISPLACEMENTS  FOR  SIMPLY  SUPPORTED  RECTANGULAR  0-90-0  PLATE  WITH  SINUSOIDAL  LOAD 


Fi:i  (:x2  Mesh)  .409  .646  .0336  .245  .286  5.166 


269  .0213  .339 


41 


TABLE  8 


Nondimensional ized  Fundamental  Frequencies  of  Simply  Supoorted, 
Square,  Cross-ply  Plate  of  Problem  1. 

(u  =  CO  -p-  V|^  ) 


Aspect 

Ratio 

Fi nite 

El ement 

Solution 

Closed  Form 
Solution 
Pcference  [4] 

QHD28 

QHD40 

warn 

2 

5.860 

5.525 

5.824 

5.500 

9,780 

9.757 

9.706 

9.359 

. 

10 

15.440 

15.340 

15.276 

15.145 

20 

17.850 

17.719 

17.628 

17.665 

25 

18.246 

18.103 

18.006 

18.093 

100 

18.964 

18.805 

18.704 

18.733 

The  Effects  of  Linear  (z).  Quadratic  (z^)  and  Cubic  (z^ 


TABLE  10 


Lamina  Stacking  Sequence 

Aspect  Ratio 
a/t 

[-45,  +45]^ 

[0,  ±45,  90], 

SELECTED  DA>UGE  CRITERIA 


TRANSVERSE  DISPLACEMENT  VS.  ASPECT 
INDRICAL  BENDING  OF  BIDIRECTIONAL  AND 
SYMMETRIC  3-PLY  LAMINATES 


NSV 


NORMALIZED  IN-PLANE  STRESS 
OF  SYMMETRIC  L 


(PRESENT) 


NORMALIZED  IN-PLANE  STRESS  FOR  CYLINDRICAL  BENDING 


PORTED 


FIGURE  11:  NORMALIZED  IN-PLANE  SHEAR  STRESS  FOR  SIMPLY 
SUPPORTED  0-90-0  SQUARE  PLATE  (x=0,  y=0,  S=4 
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FIGURE  13:  NORMALIZED  S 

SIMPLY  SUPPORTE 
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NORMALIZED  DISPLACEMENT  AND  IN-PLANE 
PECT  RATIO  (REDUCED  AND  FULL  INTEGRAT 
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DISPLACEMENTS(CM*  lOD-3) 


Z2 


30 


37 


2S 


23 


30  . 


\  7  . 


1  5  . 


1  2  . 


1  0  . 


S  . 


3  . 


0  . 


FIGURE 


65 


.5-1 


T IMC  (MiCtOSCCOMDS) 

21:  ONE-MODE  DISPLACEMENT  VS.  TIME  RESPONSE  OF  A  COARSE-MESH 

(0/90)  LAYUP  SQUARE  PLATE  UNDER  SUDDENLY-APPLIED  SINUSOIDAL 
LOADING 


DISPWCEMENTS(CM*IOD-3) 


[•lil 


TIME  (MICKOSCCONOS) 

FIGURE  22;  ONE-MODE  DISPLACEMENT  VS.  TIME  RESPONSE  OF  A  FINE-MESH 

(0/90)  LAYUP  SQUARE  PLATE  UNDER  SUDDENLY-APPLIED  SINUSOIDAL 
LOADING 


[eiKi-e], 


Failure  Curves  for  [e/9o/-e)^  Tri-directlonal  Laminates 
Due  to  In-Plane  Loading 
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APPENDIX  lA  -  HIGHER  ORDER  DISPLACEMENT  MODELS 

QHD40 

NODAL  DEGREES  OF  FREEDOM: 

Corner  Nodes  -  {Uq  Vq  Wq  %  <t>x 
Mid-side  Nodes  -  {Wq  tp^  ipy}  ^ 

DISPLACEMENT  FIELD: 

u  =  Uq  +  Zl|)x  +  z^<t>x 
V  =  Vo  +  Z4)y  +  z^<t>y 

w  =  Wo 

where; 

’Jq’  ''o’  <t'y  :  n  X  y  xy}^  {a} 

Wq.  'I'y  :  (1  X  y  x^  xy  y^  xV  xy^}^  {3} 

STRESS  FIELD: 

i.  From  constitutive  relations  -  o-j  =  Cijc-jj  (orthotropic  mat.) 

y2) 

Oyy  =  f(zS  x^,  y^) 

oxy  =  f(z^,  x^,  y^) 

ii.  From  equilibrium  considerations  -  o-ji,-:  =  0 


k* 


QHD28 


NODAL  DEGREES  OF  FREEDOM: 


''o  "^0  'I'y.  'Py  '^’x  ^y^ 


DISPLACEMENT  FIELD: 


u  +  Uq  +  ztjj^  +  z^4)x 


V  =  ZiPy  +  2^ 


W  = 


where’. 


*^0 *  ^0*  ’  ^x*  ^  ^y)  (otl 


STRESS  FIELD: 

i.  From  constitutive  relations  -  o-j  =  C-jjc-jj  (orthotropic  mat.) 


^xx 


=  f(z^  X,  y) 


Oyy  =  f(zS  X,  y) 


oxy  =  f(z^,  X,  y) 


n 


From  constitutive  considerations  -  oij.j 
oxz  =  f(z') 


=  0 


Oyz  =  f(z^) 


0^2  =  constant 


V  .*.*,■ 


•/  s 


APPENDIX  IB 


MODIFIED  KIRCHHOFF  FORMULATION 


QD32 


NODAL  DEGREES  OF  FREEDOM: 

Corner  Nodes  {wq  Vq  w  g—  Yx 

X  y 

Mid-Side  Nodes  (w 

DISPLACEMENT  FIELD: 


where; 


w  =  f(x,  y) 


/Ow  . 

U  =  Uq  -  +  -x 


y) 


Uq.  'o"  X’  y  ■  *  S' 

w  =  { 1  X  y  xy  x^y  xy‘  y^  x"  x^  xy^  y"*  xV  xy**) 


STRESS  FIELD: 


i.  From  constitutive  relations  -  (orthotropic  mat.) 

=  f(z,  x\  y^) 

Oyy  =  f(z,  x^  y^) 

oxy  =  f(z,  x^  y^) 


ii.  From  equilibrium  considerations  -  Oi  a  i  =  0 


NODAL  DEGREES  OF  FREEDOM 


‘"o  "o  «  ^-x  ’fy) 


T  ^ 


Corner  Nodes 
Mid-side  Nodes 
Center  Node 


'W  oW  Y, 


^^0  '^0  r-x  ? 
{w}-' 


{w  1^1^} 


?x  9y 


DISPLACEMENT  FIELD 
w  =  f{x,  y) 

=  "0  -  +  Yx) 

V  =  Vg  -  z  +  YxJ 

where; 

LI^,  V^,  V,  y  :  a  X  y}^( 

w  :  -.I  x  y  x'  xy  y'  x"y  xy'  y^  x"*  x^y  x’y  xy’  y} 
STRESS  FIELD 


i.  From  constitutive  relations  -  c^-  =  *'ij'  ij  (orthotropic  mat.) 

XX  '  ^  ’  y  ) 

Oyy  =  f (z,  xy  Y  ) 

o^y  =  f(z,  x’  ,  Y  ) 

i.  From  e<  lilibrium  cons  iderations  -  •iyij  =  0 


appendix  II 


HYBRID-STRESS  FORMULATION 


STRESS  FIELD: 

X  +  Bjy  +  ei^xy)  +  z(S5  +85  x  +  Sjy  +  SsXy) 

+  z-  (So  +  £5iox  +  Buy  -  •^2i^4xy) 

=  (3i2  +  3i3X  +  Bi^y  -  B.,xy)  +  z(Bi5  +  SisX  +  S^y  + 

+  z^(;3i9  -  BicX  "  BnY  p  34xy) 

“^xy  "[‘'-O  B21  +  Bii  )x  +  (-Bi  -  SiQ  “ 

+  (-B9  -  Bi9)xy]  +  z  [S23+  B24X  +  Bzsy  +  39-  ^ 

^'’■[326  B21  X  +  B22y  (pB9+  p^l^^xy] 

'^xz  ~  ~  ^  +  622  )  ■*■  Bi,y  +  (“Bo”  Bi9  )x] 

+  z2)[S(-,  +  625  +  Boy  -  -^(pg  +Bi9)x] 

+  ^(-h"-  z3)[,iio  -  p.B4y  +  B22  p(B9+  Bi9)xJ 

h  ^  Ip 

“yz  =  -  2)  [  -  J  Boi  +  -J  Bu  )  +  (-Bg  -  Big  )y  -  ;-;,x] 

1 

+  Z")l  i’2'f  -  -j^-(Bg  +  Big  )y  +  Bi7  +  B  13  X  ] 

+  ^  (-h"-  Z^)[B2i+  p^Bg  +  Bi9)y  ■  !'11+  pLoX] 

-•  =  (i:g  +  Big)  l-(h  +  z)2-  (ll^li--  +  Z^)  ,  (  Bh'lt  41^-2  ^  z 


APPENDIX  III 


MASS  MATRIX  FORMULATION 


The  mass  matrix  for  elements  under  development  is  easily  arrived  at 
by  considering  kinetic  energy  in  the  form 


=  I 


+  v^  +  w^)dV 


where  u,  v  and  w  represent  displacements,  p  is  the  mass  density  and  the 
dot  superscript  denotes  velocity.  Defining  velocities  in  terms  of  element 
shape  functions  gives 

T  =  I  j  p{Nu}[Nu]  +  (NvlCNj  +  {N„}[Nw]  dV  {A} 

V 

which  is  the  classical  form 
T  =  I  [Af  [MKaI 

The  element  mass  matrix  [M]  is,  therefore,  specified  as 

=  J  p(N^}fN^}  +  {NyliNv)  +  {N^HNwl  dV 
V 


[M] 


Note  that  the  shape  functions  [N^-]  involve  distance  from  the  mid-plane 
of  the  element  to  a  layer  denoted  by  Z  and,  therefore,  the  mass  matrix 
definition  provided  not  only  represents  mid-plane  inertial  effects  but 
also  rotatory  inertia  as  well. 


I'j 1.  I  *1/ 'J »V *  V.*  "  k* *!**-.’'  '-'-"-^'-T  *  *-* ».■  »■■  *■' 


APPENDIX  IV  -  LARGE  DISPLACEMENT  FORMULATION 


Based  on  Green's  Strain  Tensor,  the  following  procedure  is  utilized 
to  obtain  the  large  displacement  and  the  geometric  stiffness  matrices. 

Let  N  be  shape  functions  relating  displacements  at  any  point  in  the 
element  {5}  to  nodal  displacements  {A}  such  that 


{6}  =  [N]{a} 


Also  let  {N-j,j}^  denote  those  shape  functions  associated  with  the  i^^ 
displacement  field  (i  -  u,v,w)  and  ",j"  denotes  the  differentiation  with 
respect  to  the  coordinate,  i.e.,  where  Xi  =  x,  Xj  =  y  and  X3  =  z. 


Then,  the  strain  given  by 

X  X 


-  =  -^  +  i 

'XX  ^  2 


''X 


can  be  written  as 


‘XX  “  -'HPX 


2  '  !'''\l’X'’  '-'^U’X''  '-"‘V’X^  ‘''W>X"  ^'^w»x 


Siliiilarlyj  the  shear  strain  can  bo  represented  by 


c-  -!  •!  T  ^  ,T1  ..T  ,  ,.-Ti  J,.,  .  ^  J...  , 

“Xy  ["u’y'  '  ‘'^V’x-  |  '-u’X^  -’u ’y  ’ '  "v  >X'  iN^/, yi 


T 


''w’x‘  ''^W’yjj  i’’'-' 


The  strain  field  in  indicial  notation  is  expressed  by 


-  2 


;-.Ws 


•  i 


»  ..-.fl 


I _ i 


.  .  N. 

I  « 


Then  the  increinental  representation  becomes 


+  {A}"^  [{Nk.i?^f\>j^] 


•  *  ^  ^  , 


V- 


But  the  second  term  can  be  expressed  as 


{A}'’’[{Nk,j}’'’{Nk,T}]  {6A} 


Thus  combining  terms 


"  ■  "  “  ►  ' 

».  A.  i^JUL* 


{Nj,i}’'’]{6A}  +{A}1‘  [{Nk,i}’^{Nk,j}  +{Nk,j}’'’{Nk,i}]  foA} 


[Bq]  =  +  {Nj,i}'] 

[BJ  =  [t\.i)^(Nk,j}  MNk,j.}^  {Nk,i}] 


{A}'[M,,] 

{A}'[Myy] 

{A}'^[Mxy] 

{A}^[Mzz] 

{A}^[Mxz] 

{A}^[My,] 


{6A^j}  -  [Bq]{cA}  +  [BlKAA} 

where  [Bq]  is  the  linear  component  and  [Bl]  is  the  large  displacement  component 
Having  the  definitions  for  [Bq]  and  [B^],  the  small  and  large  displacement 
matrices  [Kq]  and  [K|^]  are  represented  as 

[Kq]  =  J  [Bo]^[D][BQ]dV 

V 

[Kf]  =  j  |[Bl]^[D][Bq]  +  [Bl]^[D][Bl]  +  [Bq]^[D][Bl]  [dV 


The  geometric  stiffness  matrix  is  also  derived  from  [Bl]  and  it  has  the 
following  form 


~  J" 


1  + 


Oxy[M> 


Where  the  o's  are  the  stress  components  and  again  integration  is  on  a 
layer  by  layer  basis. 


